Validation

import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
from rasterio import features

import matplotlib.pyplot as plt
import sklearn
from sklearn.metrics import confusion_matrix, f1_score, accuracy_score, classification_report

from pathlib import Path
from IPython.display import display

import plotly.express as px
import plotly.figure_factory as ff
import plotly.graph_objects as go
import plotly.offline
plotly.offline.init_notebook_mode()

print('All libraries successfully imported!')
print(f'Scikit-learn: {sklearn.__version__}')

All libraries successfully imported!
Scikit-learn: 0.24.2

Set directory

computer_path = '/export/miro/ndeffense/LBRAT2104/'
grp_letter    = 'X'

lut_path  = f'{computer_path}data/LUT/'

# Directory for all work files
work_path = f'{computer_path}GROUP_{grp_letter}/WORK/'

in_situ_path = f'{work_path}IN_SITU/'
classif_path = f'{work_path}CLASSIF/'

am_path = f'{work_path}ACCURACY_METRICS/'


Path(am_path).mkdir(parents=True, exist_ok=True)

print(f'Accuracy Metrics path is set to : {am_path}')

Accuracy Metrics path is set to : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/ACCURACY_METRICS/

Set parameters

site = 'NAMUR'
year = '2020'

feat_nb = 2


no_data = 0

reclassif_flag = True

if reclassif_flag:
    field_classif_code = 'grp_A_nb'
    field_classif_name = 'grp_A'
else:
    field_classif_code = 'grp_1_nb'
    field_classif_name = 'grp_1'

Set filenames

s4s_lut_xlsx   = f'{lut_path}crop_dictionary_new.xlsx'

in_situ_val_shp = f'{in_situ_path}{site}_{year}_IN_SITU_ROI_VAL.shp'
in_situ_val_tif = f'{in_situ_path}{site}_{year}_IN_SITU_ROI_VAL.tif'

if reclassif_flag:
    classif_tif = f'{classif_path}{site}_{year}_classif_RF_feat_{feat_nb}_{field_classif_name}.tif'
else:
    classif_tif = f'{classif_path}{site}_{year}_classif_RF_feat_{feat_nb}.tif'


cm_csv = f'{am_path}{site}_{year}_CM.csv'
cm_html = f'/export/miro/ndeffense/LBRAT2104/GIT/eo-toolbox/notebooks/7_Classification/figures/{site}_{year}_CM.html'
am_html = f'/export/miro/ndeffense/LBRAT2104/GIT/eo-toolbox/notebooks/7_Classification/figures/{site}_{year}_AM.html'

print(f'Classification file used : {classif_tif}')
print(f'Validation polygons used : {in_situ_val_shp}')

Classification file used : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/CLASSIF/NAMUR_2020_classif_RF_feat_2_grp_A.tif
Validation polygons used : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/IN_SITU/NAMUR_2020_IN_SITU_ROI_VAL.shp

1. Prepare data

1.1 Rasterize in situ data validation shapefile

print(f'Raster template file : {classif_tif}')

# Open the shapefile with GeoPandas

in_situ_gdf = gpd.read_file(in_situ_val_shp)

# Open the raster file you want to use as a template for rasterize

src = rasterio.open(classif_tif, "r")

# Update metadata

out_meta = src.meta
out_meta.update(nodata=no_data)

crs_shp = str(in_situ_gdf.crs).split(":",1)[1]
crs_tif = str(src.crs).split(":",1)[1]

print(f'The CRS of in situ data is    : {crs_shp}')
print(f'The CRS of raster template is : {crs_tif}')

if crs_shp == crs_tif:
    print("CRS are the same")

    print(f'Rasterize starts : {in_situ_val_shp}')

    # Burn the features into the raster and write it out

    dst = rasterio.open(in_situ_val_tif, 'w+', **out_meta)
    dst_arr = dst.read(1)

    # This is where we create a generator of geom, value pairs to use in rasterizing

    geom_col = in_situ_gdf.geometry
    code_col = in_situ_gdf[field_classif_code].astype(int)

    shapes = ((geom,value) for geom, value in zip(geom_col, code_col))

    in_situ_arr = features.rasterize(shapes=shapes,
                                     fill=no_data,
                                     out=dst_arr,
                                     transform=dst.transform)
    
    dst.write_band(1, in_situ_arr)

    print(f'Rasterize is done : {in_situ_val_tif}')

    # Close rasterio objects
    src.close()
    dst.close()

else:
    print('CRS are different --> repoject in-situ data shapefile with "to_crs"')

Raster template file : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/CLASSIF/NAMUR_2020_classif_RF_feat_2_grp_A.tif
The CRS of in situ data is    : 32631
The CRS of raster template is : 32631
CRS are the same
Rasterize starts : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/IN_SITU/NAMUR_2020_IN_SITU_ROI_VAL.shp
Rasterize is done : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/IN_SITU/NAMUR_2020_IN_SITU_ROI_VAL.tif

1.2 Build y_pred and y_true

# Open in-situ used for validation
src = rasterio.open(in_situ_val_tif, "r")
val_arr = src.read(1)
src.close()

# Open classification map
src = rasterio.open(classif_tif, "r")
classif_arr = src.read(1)
src.close()

# Get the postion of validation pixels
idx = np.where(val_arr == no_data, 0, 1).astype(bool)

# Ground truth (correct) target values
y_true = val_arr[idx]

# Estimated targets as returned by a classifier.
y_pred = classif_arr[idx]

Sometimes, some classes do not appear in the classification map, they are not predicted by the Random Forest.

This means that some classes in y_true don't appear in y_pred.

classes_all  = sorted(np.unique(y_true))
classes_pred = sorted(np.unique(y_pred))

classes_missing = set(y_true) - set(y_pred)

print(f'{len(classes_missing)} classes are missing in the classification (y_pred) : {classes_missing} \n')

print(f'All training classes :\n {classes_all}')
print(f'All predicted classes (at least once) :\n {classes_pred}')

0 classes are missing in the classification (y_pred) : set() 

All training classes :
 [3, 6, 8, 9, 17, 21, 22, 111, 112, 115, 117, 119, 121, 143, 151, 181, 192]
All predicted classes (at least once) :
 [3, 6, 8, 9, 17, 21, 22, 111, 112, 115, 117, 119, 121, 143, 151, 181, 192]

1.3 Find class names

lut_df = pd.read_excel(s4s_lut_xlsx)

classes_name = lut_df[lut_df[field_classif_code].isin(classes_all)].sort_values(field_classif_code)[field_classif_name].drop_duplicates().to_list()

for code,name in zip(classes_all, classes_name):
    print(f'{code} - {name}')

3 - Grassland and meadows
6 - Forest
8 - Build-up surface
9 - Water bodies
17 - Leguminous crops
21 - Fruits trees
22 - Vineyards
111 - Wheat
112 - Maize
115 - Barley
117 - Oats
119 - Other cereals
121 - Leafy or stem vegetables
143 - Other oilseed crops
151 - Potatoes
181 - Sugar beet
192 - Fibre crops

2. Confusion Matrix

2.1 Compute Confusion Matrix

cm = confusion_matrix(y_true, y_pred)

cm_df = pd.DataFrame(cm)

cm_values = cm_df.values

cm_df.columns = classes_all
cm_df.index = classes_all

cm_df.to_csv(cm_csv, index=True, sep=',')

display(cm_df)

	3	6	8	9	17	21	22	111	112	115	117	119	121	143	151	181	192
3	17445	298	20	0	5	56	20	62	0	1	3	8	0	0	0	0	0
6	37	746	30	0	0	0	0	0	0	0	0	0	0	0	0	0	0
8	52	14	422	13	0	3	0	20	0	2	0	0	0	0	0	0	0
9	0	0	26	30	0	0	0	0	0	0	0	0	0	0	0	0	0
17	151	0	7	0	1116	0	0	20	0	0	0	0	0	0	1	0	8
21	686	0	0	0	2	0	0	0	0	0	0	0	0	0	0	0	0
22	14	0	0	0	0	40	0	0	0	0	0	0	0	0	0	0	0
111	90	3	47	0	14	0	4	20100	0	1059	65	750	0	2	0	1	95
112	28	0	25	0	2	0	0	1	6344	0	3	10	4	7	87	22	2
115	199	0	4	0	0	0	0	1509	0	2726	136	401	0	21	0	0	0
117	1	0	0	0	0	0	0	11	0	0	88	0	0	0	0	0	142
119	12	0	11	0	0	0	0	4005	0	0	112	455	0	0	0	0	45
121	5	0	0	0	0	0	0	0	738	0	0	0	23	0	0	19	0
143	51	0	0	0	0	0	0	18	0	24	0	0	0	557	0	0	0
151	0	1	6	0	0	0	1	0	271	0	0	4	6	0	3538	1035	0
181	2	14	1	0	4	2	1	0	389	0	0	0	166	0	21	4950	0
192	0	0	2	0	0	0	0	605	0	0	40	45	0	0	0	0	59

2.2 Plot Confusion Matrix

# invert z idx values
z = cm[::-1]

x = classes_name
y = x[::-1].copy() # invert idx values of x

# change each element of z to type string for annotations
z_text = [[str(y) for y in x] for x in z]

# set up figure 
fig = ff.create_annotated_heatmap(z,
                                  x=x,
                                  y=y,
                                  annotation_text=z_text,
                                  colorscale='spectral',
                                  reversescale=True)

# add title
fig.update_layout(title_text=f"Confusion Matrix - {site}, {year}")

# adjust margins to make room for yaxis title
#fig.update_layout(margin=dict(t=200, l=200))

#fig.update_xaxes(tickfont_size=20)
#fig.update_yaxes(tickfont_size=20)
#fig.update_layout(font_size=25)

# add colorbar
#fig['data'][0]['showscale'] = True

fig.show()
fig.write_html(cm_html, full_html=False)

3. Accuracy metrics

3.1 Compute Accuracy Metrics

If you decide that you are not interested in the scores of classes that were not predicted, then you can explicitly specify the classes you are interested in (which are labels that were predicted at least once).

acc_metrics_str  = classification_report(y_true,
                                         y_pred,
                                         target_names=classes_name,
                                         labels=classes_all,
                                         digits=3)

print(acc_metrics_str)

                          precision    recall  f1-score   support

   Grassland and meadows      0.929     0.974     0.951     17918
                  Forest      0.693     0.918     0.790       813
        Build-up surface      0.702     0.802     0.749       526
            Water bodies      0.698     0.536     0.606        56
        Leguminous crops      0.976     0.856     0.913      1303
            Fruits trees      0.000     0.000     0.000       688
               Vineyards      0.000     0.000     0.000        54
                   Wheat      0.763     0.904     0.827     22230
                   Maize      0.819     0.971     0.889      6535
                  Barley      0.715     0.546     0.619      4996
                    Oats      0.197     0.364     0.255       242
           Other cereals      0.272     0.098     0.144      4640
Leafy or stem vegetables      0.116     0.029     0.047       785
     Other oilseed crops      0.949     0.857     0.901       650
                Potatoes      0.970     0.728     0.832      4862
              Sugar beet      0.821     0.892     0.855      5550
             Fibre crops      0.168     0.079     0.107       751

                accuracy                          0.807     72599
               macro avg      0.576     0.562     0.558     72599
            weighted avg      0.774     0.807     0.783     72599

3.2 Compute Overall Accuracy

oa = accuracy_score(y_true, y_pred)
oa = round(oa*100, 2)

print(f'Overall Accuracy : {oa}%')

Overall Accuracy : 80.72%

3.3 Plot Accuracy Metrics

acc_metrics_dict = classification_report(y_true, y_pred,target_names=classes_name, output_dict=True)

am_df = pd.DataFrame.from_dict(acc_metrics_dict).round(3)

am_df = am_df.iloc[:,:-3]

#am_df = pd.concat([am_df, nb_df])

#am_df = am_df.sort_values(by='pix_count', ascending=False, axis=1)

class_name = am_df.columns.to_list()
precision = am_df.loc['precision'].to_list()
recall    = am_df.loc['recall'].to_list()
f1_score  = am_df.loc['f1-score'].to_list()

fig = go.Figure(data=[
    go.Bar(name='Precision', x=class_name, y=precision, text=precision, textposition='auto'),
    go.Bar(name='Recall', x=class_name, y=recall, text=recall, textposition='auto'),
    go.Bar(name='F1-score', x=class_name, y=f1_score, text=f1_score, textposition='auto')
])
# Change the bar mode
fig.update_layout(title_text=f'Accuracy Metrics - {site}, {year}',
                  barmode='group')

#fig.update_xaxes(tickfont_size=30)
fig.update_yaxes(tickfont_size=10, range=[0,1])
#fig.update_layout(xaxis_title=None, font_size=10)
#fig.update_layout(legend=dict(font=dict(size=25)))

fig.show()
fig.write_html(am_html, full_html=False)